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Abstract. 

We present an algorithm to evaluate large deviation functions associated to history- 
dependent observables. Instead of relying on a time discretisation procedure to 
approximate the dynamics, we provide a direct continuous-time algorithm valuable 
for systems with multiple time scales, thus extending the work of Giardina, Kurchan 
and Peliti [1]. 

The procedure is supplemented with a thermodynamic-integration scheme which 
improves its efficiency We also show how the method can be used to probe large 
deviation functions in systems with a dynamical phase transition - revealed in our 
context through the appearance of a non-analyticity in the large deviation functions. 
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1. Introduction 

The statistical physics of equilibrium systems was first designed to reproduce the 
macroscopic predictions of thermodynamics, but it was soon realised [2] that it also 
provides a well-suited frame to describe the fluctuations of physical observables. Such 
a theory is not available when dealing with non-equilibrium systems or dynamical 
observables, as one lacks thermodynamics functions such as the free-energy. Over the 
last decade, there has been a growing interest within the physics community in the theory 
of large deviation functions [3], as it appeared that they could fill this gap in some cases. 
For instance, the fluctuations of particle or energy currents Q flowing through a system 
in the long time limit can be obtained from the associated large deviation function 
7r(q = Q/t): 

Prob(g, t) ~ e t7r(q) as t -> oo (1) 

The function n(q) is a dynamical analog of the intensive entropy in the microcanonical 
ensemble and the long time limit plays the role of the thermodynamic limit. In the 
past few years, a huge amount of effort has been devoted to the study of a strikingly 
simple symmetry of the large deviation function of the injected power, the so-called 
fluctuation theorem. The results obtained range from theoretical and numerical studies 
[4, 5, 6, 7, 8, 9] to experimental applications (see [10] for a glimpse at the literature). 

Beyond its symmetries, ir(q) itself expectedly bears information on the current 
flowing through the system. For instance, the large deviation function n(q) can be fully 
determined in diffusive systems [11, 12] and its scaling properties differ from those of 
superdiffusive ones (see [13] for explicit results). In the context of dynamical systems 
also, large deviation functions associated with more general observables have been 
introduced [14] and their determinations has received a broad interest [16, 18, 19, 17]. 
Although the variety of results obtained from these approaches raises hopes for an 
out-of-equilibrium thermodynamics, the determination of large deviation functions is in 
general a hard task to achieve and most exact results are confined to simple systems 
or peculiar models [11, 12, 13]. In more complex cases, we have to rely on numerical 
evaluations. 

However, the very definition of large deviations renders their direct numerical 
observation almost impossible, as the probability to observe a value of Q/t far from its 
average typically decreases exponentially with time. Giardina, Kurchan and Peliti [1] 
introduced a numerical procedure which overcomes this difficulty for discrete time 
Markov chains. Nevertheless, most physical processes evolve continuously in time, and 
one thus has to choose an arbitrary time step dt to discretise the dynamics, balancing 
between algorithm efficiency and errors arising from the approximation. Typically, dt 
has to be smaller than any time scale of the system, yet too small a value only increases 
the simulation duration, since most of the time steps would then be spent in rejected 
moves. This issue already affects the standard Monte Carlo algorithms, and becomes 
quite unavoidable when dealing with the computation of large deviation functions. 
Indeed, even systems featuring a single time scale in the steady state present in general 
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different time scales in their large deviations, depending on the kind of histories probed, 
which makes the choice of dt strenuous. A simple example is given by traffic flow models, 
in which the typical time scale is fixed on average, but vary by a factor equal to the 
number N of "cars" when comparing jammed histories (where 0(1) cars move) and free 
flowing histories (where O(N) cars move). 

In the present paper we propose a new procedure where the time discretisation 
issue is bypassed with a direct continuous-time approach. The outline of the paper 
is as follows: in section 2, we recall the continuous-time formalism, define the large 
deviation functions and present the algorithm. In section 3, we study three systems 
where the continuous-time approach proves useful: the symmetric simple exclusion 
process, its asymmetric, out of equilibrium counterpart, and the contact process, for 
which a dynamical phase transition occurs. 



2. Formalism and Algorithm 

2.1. Continuous time Markov chains 

We consider a system described by a finite number of configurations {C}, whose evolution 
is determined by the transition rates W(C — > C) between different configurations. The 
probability P(C,t) to find the system in C at time t evolves according to the master 
equation 

d t P(C, t) = J2 W(C -> C)P(C, t) - r(C)P(C, t) (2) 

where the escape rate r(C) reads 

r(C) = £V(C-C'). (3) 

Starting from C , the system evolves through a succession of configurations {Ck}o<k<K, 
jumping from C fc to Ck+i at time t k with probability W ^^ k+1 ^ (See Figure 1). Note 
that contrary to the discrete time case, Ck and Ck+i are necessarily different. The time 
elapsed between two consecutive jumps is a random variable: if the system arrives at 
time t in the configuration C , the next move occurs at a time t±, distributed according 
to a Poisson law: 

p(t 1 \C ,t )=r(C )e-^' t ^ (4) 



2.2. Large deviation functions 

Let us consider an observable A extensive in time, which can be written as a sum along 
a history {Ck}o<k<K of elementary contributions ac k c k+1 ' 

A= a c k c k+l (5) 

0<fc<A"-l 

This form is quite generic and most of the commonly studied observables fall in this 
class. For instance, if A is the overall current of a one dimensional lattice gas, ace 
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Figure 1. An history of the system between time to = and time t. 

is the contribution of a single particle jump (See Section 3.1). Moreover, to compute 
the average of a static observable O along a history, one simply takes ace — 0(C) an d 
recover (O) = j. 

In equilibrium statistical mechanics, the difficult task of computing the entropy 
can be conveniently substituted by the determination of the free energy. Here, instead 
of working at fixed value of A (microcanonical ensemble) to compute n(A/t), one 
rather introduces a parameter s which fixes the average of A (canonical ensemble), s is 
intensive in time and plays a role analogous to the inverse temperature in equilibrium 
thermodynamics. This leads us to introduce the dynamical partition function 

Z(s,t) = (e' sta ) ~ e'^ (s) as t -> 00 (6) 

where the average (. . .) is taken over all histories between and t and a — A/t is intensive 
in time. The large deviation function iPa(s) of Z(s,t) is the Legendre transform of 11(a) 

4>a(s) = max [71(a) — sa] . (7) 

a 

Under quite general conditions [3], this relation can be inverted and one can get 11(a) 
knowing iPa(s) through: 

n(a) = max [iPa(s) + sa] (8) 

s 

To compute Z(s, t), let us first write the master equation obeyed by the joint probability 
P(C, A, t) of being in configuration C with a value A for the observable at time t: 

d t P(C, A,t) = J2 W ( C ' C ) p ( c> > A ~ a cc, t) - r(C)P(C, A, t) (9) 

The Laplace transform P(C,s,t) = ^ A e~ sA P(C, A,t) then evolves with 

d t P(C, s,t) = J2 W s(C - C)P(C, s, t) - r(C)P(C, s, t) (10) 

C'^C 

where the s-modified rates W s are given by 

W S (C ->C) = e- sa ccW(C ->C) (11) 

From equation (10), one sees that Z(s,t) = P(C, s,t) behaves at large time as 
e <pA(s)t yjr^Qj-Q ip A (s) is the largest eigenvalue of a (not probability conserving) evolution 
operator, which justifies the asymptotic behaviour in (6). The determination of the 
large deviation functions iPa(s) then amounts to the computation of this eigenvalue, 
which we address in the next section. 
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2.3. A cloning algorithm 

Let us consider the s -modified Markov dynamics defined by the rates W s , whose evolution 
operator reads 

(W s ) cc , = W S (C -> C) - r s (C)5 cc > (12) 

where 

r s (C) = J2w s (C^C). (13) 
c 

The evolution of P(C, s, t) can be written as 

d t P(C, s, t) = ( W *W p (C, s, t) + [r s (C) - r(C)]P(C, s, t) (14) 
c 

The corresponding dynamics alternates changes of configuration determined by the s- 
modified rates and exponential evolution of the "non-conserved probability" P(C, s, t) 
with rate r s (C) —r(C) (corresponding respectively to the first and second term of the r.h.s 
of (14)). More details are given in Appendix A. The evolution of Z(s, t) = J2 C P(C, s, t) 
is consequently represented by a population dynamics a la Diffusion Monte Carlo [20]. 
Let us consider Mo clones of the system evolving in parallel with the s modified dynamics. 
At any change of configuration of any clone c a , 

(1) c a jumps from its configuration C to another configuration C with probability 
W s (C^C')/r s (C). 

(2) The time interval At until the next jump of c a is chosen from the Poisson law (4) 
of parameter r s (C). 

(3) The clone c a is either cloned or pruned with a rate y(C) = e A *^ 3 ^ c '^ _r ^ c '^ 

a) One computes y = \_y(C) + e\ where e is uniformly distributed on [0, 1]. 

b) If y = 0, the copy c a is erased. 

c) If y > 1, we make y — 1 new copies of c a 

Such a cloning procedure modifies the total number of clones by a factor X = J ^ + J ^~ 1 , 
which represents the exponential evolution of P(C, s,t). Finally, Z(s,t) is simply given 
by the increase of the population: 

Z(s,t) = ^ (15) 

However, such an algorithm may result in an exponential growth or a complete decay 
of the clones and we consequently add a 4 th step to maintain their number constant: 

(4) If y = 0, one clone eg ^ c a is chosen at random and copied, while if y > 1, y — 1 
clones are chosen uniformly among the N + y—l clones and erased. To reconstruct 
Z(s,t), we keep track of all X factors 

The large deviation function iPa(s) is then recovered from the long time behaviour of 
the product of the cloning factors: 

-\nX 1 ...X K = -ln(e- sta ) ~ ^ A (s) as t -> oo (16) 
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where /C is the total number of configuration changes among all the histories between 
and t. 

Another interpretation of such a population dynamics has been discussed by 
Grassberger in [15] for the discrete time case and we shall present it shortly. A formal 
integration of (14) leads to 

(e~ Sta ) = (^ e Io r ^(t))-r(C(t))dt^ ( 17 ) 

where the average (. . .) s is taken over all trajectories from to t of the s-modified 
dynamics. The change of rates W — > W s can be seen a sequential importance sampling, 
which favours histories relevant for the computation of Z(s,t). The exponential within 
the r.h.s. of (17) can be seen as a weight over trajectories. Along the simulation, clones 
with a high weight wh are replaced by wh clones with a weight 1, while clones with low 
weight wl are either deleted with probability 1—wl or given a weight 1 with probability 
wl- This is very close to the cloning strategy followed at step (3). 

Let us finally note that one can define a new measure over the space of trajectories, 
such that the average of an observable B reads 

W^^l (18) 

This is the measure observed in the simulations with a fixed number of clones, as one 
can see that the denominator is nothing but the unconstrained population at time t. 

2.4- Thermodynamic integration 

The direct computation of iPa(s) is in general quite noisy, because of the finiteness of the 
clones population. One can alternatively compute ip' A (s) and then integrate it. From 
the definition of ipA '■ 

Ms) = ]^(e- Sta ) (19) 
one gets by derivation 

= J jf^ = (20) 
which is nothing but the average value of a within the population of clones. Finally, 

Ms) = - [ ^dr (21) 
Jo 

Thanks to the integration, the noise is smoothed out. 
3. Three examples 

3.1. The Symmetric Exclusion Process (SEP) 

We now apply the algorithm to the large deviations of the total current Q in the 
Symmetric Exclusion Process [21] with closed boundary conditions. The system is 
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composed of iV particles diffusing on a one-dimensional lattice of size L. Each particle 
can jump with rate 1 to any neighbouring site, provided it is empty. The total current Q 
increases or decreases by 1 at every move, depending on the direction of the jump. 
Using the notation (5), ace = 1 or —1 when a particle moves to its right or to its left, 
respectively. 
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Figure 2. Numerical evaluation of ^q{s) for the Simple Exclusion Process (N = 200, 
L = 400). (a: Left) Comparison between direct numerical measurement (blue crosses) 
and result from thermodynamic integration (red circles). 

(b: Right) Comparison between numerical results (red circles) and analytical 
prediction (22) valid for small s (blue line). 

As in other exclusion processes [9, 13], the large deviation function ?Pq(s) is 
extensive in the system size and we rather study the rescaled large deviation function 
in the thermodynamic limit (N and L are large, the density p = N/L being 
fixed). Though on average the total current is zero (the moves are symmetric), the 
variance is finite and iPq(s) reads for small s [11, 21] 

^ Q (s)= p (l- p )s 2 + 0(L S *) (22) 

In this regime, the fluctuations are Gaussian and our algorithm yields results in perfect 
agreement with the expansion (22) (see Figure 2b). For larger values of s, we lack 
an analytic expression for i/jq(s) but the algorithm still applies successfully. We found 
non-Gaussian fluctuations (Figure 2a), which correspond to very large deviations of the 
current. 

The direct measurement of iPq{s) is also compared with the results obtained from 
thermodynamic integration (Figure 2a). Both evaluations coincide within numerical 
accuracy, though the latter requires duration more than one order of magnitude smaller 
to converge. 

3.2. The Asymmetric Simple Exclusion Process (ASEP) 

We now consider the large deviations of the total current Q in a non-equilibrium system, 
the Asymmetric Simple Exclusion Process [21] with closed boundary conditions. The 
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Figure 3. Plot of the large deviation function -tiPq(s) of the Asymmetric Simple 
Exclusion Process, for L = 400 sites and N = 200 particles. The jump rates are 
p = 1.2 and q — 0.8, whence E ~ —0.2. Blue crosses and red circles correspond 
to direct computation and thermodynamic integration respectively. The asymmetry 
appears when comparing the extreme points s = ±9.5. 




25 50 75 100 25 50 75 100 
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Figure 4. (a: Left) Average profile p for s = 0.3. To minimise the overall current, the 
system develops an asymmetric shock, where only the front particles can jump easily, 
(b: Right) A typical configuration for \s\ S> E. The particles are distributed almost 
uniformly. 



system is composed of N particles diffusing on a one-dimensional lattice of size L. Each 
particle can jump with rate p to its left and q to its right, provided the arrival site 
is empty. The total current Q is defined as in the previous section. For p ^ q, a 
non-zero steady current flows through the system. We report in Figure 3 the large 
deviation function ^q(s). One notes the symmetry around E — |ln| , guaranteed by 
the fluctuation theorem. The histories which contribute to iPq{s) around s ~ E (Q « 0) 
clearly display shocks (Figure 4a) while a large current (|s| ^> E) is provided by uniform 
profiles (Figure 4b). 
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Figure 5. (a: Left) Plot of the large deviation function j^^k{s) associated to the 
number of events K in the Contact Process in a field (L — 120 sites) (b: Right) 
The dynamical phase transition occurs at s c ~ 0.057. This is exemplified by plotting 
ip' K (s) — \{K) S for different system sizes (L — 4 in black, 8 in red, 15 in blue and 50 
in magenta) 

3.3. The Contact Process (CP) 

We now turn our attention to the Contact Process [23] in one dimension. The model is 
defined on a lattice of L sites with periodic boundary conditions. Each site % is either 
empty (n$ = 0) or occupied by a particle (rt, = 1). The dynamics is defined as follows: 
particles annihilate with rate 1, while empty sites i get occupied with rate 

Willi = -> ni = 1) = A(wi_i + n i+1 ) + h (23) 

where A and h are positive constants. Note the presence of a spontaneous rate of 
creation h. When h = 0, the systems always reaches an absorbing empty state in finite 
size [24, 25], while the steady state is active for A > 1 in the thermodynamic limit. 

On the contrary, the addition of a small field h leads to an active equilibrium state 
for all A. In the mean field version, the equilibrium dynamics is still influenced by the 
presence of an inactive state. This can be seen through the study of the large deviations 
of the number of events K, a quantity which simply counts the number of configuration 
changes during an history of the system. It was proved that jr?p K (s) is non-analytic 
at a critical value s c , which goes to zero with h. As pointed out in section 2.2 the 
large deviation function iPk{s) plays the role of a dynamical free energy, whose non- 
analyticities are synonymous of dynamical phase transitions. In physical terms, this 
corresponds to the existence of two distinct classes of histories [19]: the active phase 
dominates the steady state, while large deviations corresponding to s > s c are dominated 
by inactive histories. For s = s c , the two phases coexist in a first-order fashion. The 
question whether this transition is still present in finite dimension is still open. 

Using our algorithm, we obtain evidence that this is indeed the case in dimension 
1. In Figure 5a, we plot the function iPk{s) for a large system size and for values of 
the parameters A = 3.5, h = 0.1. The two branches of the function correspond to 
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the two dynamical phases mentioned above. As in the mean field version, the non- 
analycity appears as a jump in the first derivative of ^k(s) in the thermodynamic 
limit. Using the relation i/j'k(s) = — \K S (see section 2.4), we can study the finite size 
scaling of j;V4:( s ) (Figure 5). The results support the presence of a phase transition 
at s c ~ 0.057. Around s c , the dynamics presents a superposition of "more active" and 
"less active" histories, for which the continuous-time approach is particularly helpful. 



4. Conclusions 



We have presented a simple algorithm to evaluate large deviation functions in 
continuous-time Markov chains without relying on any time discretisation. We have 
shown on specific examples that the method can be used successfully in systems where 
the presence of different time scales renders the discrete-time approach difficult. 

In the context of quantum simulations, an approach to continuous-time versions of 
Diffusion Monte Carlo was proposed by Syljuasen [26], using a continuous formalism but 
a discrete time implementation. One can expect that the algorithm we have introduced 
could also be interesting in such context. 
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Appendix A. 

An explicit formula for Z(s,t) can be written as the sum over C of the solution of (14): 

^ t) = e e f d ^ K r Wtx-! ■■■ r d^e-^-w (A.i) 

K>0C U -,Ck-i,C 

where 

r s (C) = E W S (C - C) (A.2) 

is the escape rate from the configuration C in the s-modified dynamics, 

df4 k = dt k r s (C k ) exp[-(t k - t k ^i)r s (C k )} (A. 3) 
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represents the probability of the time intervals between jumps, and 

y(Ck) tk+1 ~ tk = e^ fc+1 ~ ifc ^ r ^ Cfc ^~ r( ^ Cfc ^ (A. 4) 

is the exponential increase of P between t fc and t k+ i. We can read in this formula the 
direct expression of the cloning factors used at the step (3) of the algorithm. 
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